Dynamical properties of Au from tight-binding molecular-dynamics simulations 
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We studied the dynamical properties of Au using our previously developed tight-binding method. 
Phonon-dispersion and density-of-states curves at T=0 K were determined by computing the 
dynamical-matrix using a supercell approach. In addition, we performed molecular-dynamics simu- 
lations at various temperatures to obtain the temperature dependence of the lattice constant and of 
the atomic mean-square-displacement, as well as the phonon density-of-states and phonon-dispersion 
curves at finite temperature. We further tested the transferability of the model to different atomic 
environments by simulating liquid gold. Whenever possible we compared these results to experi- 
mental values. 



I. INTRODUCTION 

Over the last two decades atomistic simulations have 
become an increasingly important tool for modeling in 
many areas of condensed-matter physics and material sci- 
ence. The most challenging problem in computer-based 
nano-scale simulations of real materials is to find an accu- 
rate and transferable model for the atomic interactions 
that reproduces the energetic and electronic properties 
of the material. A whole hierarchy of models for atomic 
interactions have been developed, ranging from simple 
empirical potentials to sophisticated first-principles cal- 
culations based on density- functional-theory (DFT). Al- 
though DFT methods are very accurate and have been 
successfully applied to the study of a broad range of ma- 
terials and systems, they are computationally very de- 
manding. Even with today's state-of-the-art computers, 
DFT simulations with more than 100 atoms are chal- 
lenging. Empirical potentials, on the other hand, are 
less demanding and have been used to simulate systems 
with millions of atoms. This advantage is however to be 
weighed against a loss in accuracy and transferability. 

Several empirical potential methods have been used in 
the past to simulate metallic systems: the embedded- 
atom method, the effective-medium theory, Finnis- 
Sinclair potentials and the second-moment approxima- 
tion to the tight-binding modelE The decade has seen the 
emergence of a method that lies between first-principles 
and empirical potentials: the so-called tight-binding 
(TB) molecular-dynamics method. It is more accurate 
than the empirical potential methods because it explic- 
itly describes the electronic-structure of the system. TB 
is roughly three orders of magnitude faster than DFT 
based methods due to the much smaller size of the sec- 
ular equation, which makes the N"^ issue more tolerable. 
The TB njupthod has been used to study a broad range of 
materials. □ 

Recently the NRL group proposed an alternative for- 
mulation of the TB mcthOjd, which was shown to work 
well for transition metals,Ll simple metalsjj and semi- 



conductors. El This approach has been successful in deter- 
mining static properties such as structural energy differ- 
ences, elastic constants, vacancy formation energies and 
surface energies. 

Although static calculations are very useful for deter- 
mining many fundamental properties of materials, such 
calculations are limited to properties at T=0 K. Most 
problems in real materials involve processes that occur 
at finite temperature. The purpose of the present work 
is to demonstrate that our TB model can successfully be 
applied to the study of the dynamical and finite temper- 
ature properties of a representative-material, gold. Our 
previous TB parmetrization of golcH was highly success- 
ful in predicting structural properties. We have improved 
upon this parametrization in this paper. This provides 
us with an ideal test case for demonstrating the power of 
the method. 

We tested our TB parameters by calculating the elas- 
tic properties and comparing to first-principles calcula- 
tions and experiment. We also found the the phonon- 
dispersion curves and density-of-states (DOS) at T=0 K 
by calculating the dynamical-matrix using a supercell 
method.cl In addition, we performed molecular-dynamics 
(MD) simulations at various temperatures to obtain the 
temperature dependence of the lattice constant and of 
the atomic mean-square-displacement, as well as the 
electronic and phonon DOS and the phonon-dispersion 
curves at finite temperatures and a simulation of the liq- 
uid phase. Whenever possible we compare these results 
to experimental data. 

II. TECHNICAL DETAILS 

Fitting procedure for Au 

Details about our TB model can be found in Ref.i. 
In this paper we used a new TB parametrization for 
Au,lI which works well even at very small interatomic 
distances. The parameters of the model are fitted to re- 
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produce data from DFT calculations: band structures 
and total energy as a function of volume for fee, bee 
and simple cubic (sc) structures. In the present case 
the database included ten (10) fee structures, six (6) bee 
structures, and five (5) se structures. The calculations 
were performed using the general potential-.Linearized 
Augmented PrUne Wave (LAPW) method^u using the 
Perdew-Wangli-jparametrization of the Local Density 
Approximation.O In addition care was taken to include 
energies at very small volumes (down to 60% of the equi- 
librium volume) in the fitting database. This turned out 
to be very important in order to have parameters that 
could be used in the wide range of interatomic distances 
that occur during MD simulations. Finally it should be 
stressed that no experimental data is used to determine 
the parameters of the model. 

DoD-TBMD eode 

Except as noted, the results presented in this pa- 
per were obtained using the DoD-Parallel Tight-Binding 
Molecular-Dynamics (TBMD) code developed as part 
of the Computational Chemistry and Materials Science 
(CCM) contribution to the Common HFC Software Sup- 
port Initiative (CHSSI). This program was written with 
the goal of performing molecular-dynamics simulations 
of metallic systems. Althangh initially written to run 
with our TB-JIamiltonian,El this code is in fact model 
independent .113 The electronic structure is calculated us- 
ing either a 0(N'^) method such as diagonalization, or 
by using an O(N^) mfjthod called the Kernel Polyno- 
mial Method (KPM). yy The code has been written for 
both scalar and parallel computers. The parallel parts 
of the code have been written using a message-passing 
programming model r|eijying on the MPI library to deal 
with communications 

Simulation details 

To compute the dynamical-matrix we used an fee su- 
percell of 1331 atoms, obtained by replicating a primitive 
fee cell 11 times along the three primitive lattice vec- 
tors. Periodic boundary conditions are applied through- 
out this work. In the MD simulations the system consists 
of an fee supercell of 343 atoms, obtained by replicat- 
ing a primitive fee cell 7 times along the three primitive 
lattice vectors. The Brillouin-zone (BZ) is sampled us- 
ing the F-point. We checked that this is a reasonable 
approximation even for a metal: the lattice constant, 
bulk modulus and elastic constants obtained from the 
343 atom supercell and F-point sampling are within 10% 
of the values obtained using a primitive cell and a well 
converged k-point set. The MD simulations were started 
with atoms arranged on an fee lattice and random ve- 
locities drawn from a Boltzmann distribution for a tem- 
perature 2T. The MD simulation was performed in the 
micro-canonical ensemble, so at equilibrium the temper- 
ature of the lattice averaged to T, and the "potential en- 
ergy" of the system was raised by an amount 3/2 NksT, 



where N is the number of atoms and ks is Boltzmann's 
constant. The equations of motion were integrated us- 
ing the Verlet algorithm and a time step of A< = 2 fs, 
giving a total-energy conservation within IS.E / E = 10^^. 
The system was equilibrated at the desired temperature 
for 1500 time steps (3 ps). Typically another 1500 addi- 
tional steps were performed to calculate time averages. 
The finite size of the simulation cell and the finite num- 
ber of time steps in the simulation implies that the in- 
stantaneous temperature of the cell, computed from the 
kinetic energy, fluctuated around some average tempera- 
ture, which was not necessarily the target temperature. 
For the simulations conducted at a target temperature 
of 300K we found that the average temperature after the 
system reached equilibrium was 30 IK, with a standard 
deviation of 8K. At 600K, we found the average temper- 
ature to be 604K with standard deviation 20K, and at 
1200K we found 1212K and 113K, respectively. 

For the computation of the temperature dependence of 
the lattice constant we used a 64 atom fee supercell and 
sampled the BZ with four fc-points. For each volume aiid 
temperature we ran a Langevin dynamics simulationt£l 
for 2.5 ps, using a time step of 5 fs and a friction param- 
eter 7=0.05 fs-i. 

III. RESULTS AND DISCUSSION 

Equation of State 

In Fig we present the energy versus volume curves 
for a selection ofrejystal structures, calculated using the 
static TB code.EI The LAPW total energies for the fee, 
bee, and se structures, which were used in the fit, are 
also plotted on the graphs. We find that the TB Hamil- 
tonian predicts that the equilibrium fee structure has a 
lower energy than any other structure yet tested, con- 
sistent with experiment and first-principles calculations, 
and confirming the robustness of the Hamiltonian. 

Elastie properties 

The elastic properties of bulk fee Au were calculated 
using our TB parameters by means of the standard finite 
strain methodllSo and the static code. The results are 
shown, in Table I, along with comparisons to experimental 
datao and the results of first-principles LAPW calcula- 
tions. The latter calculations were also pjejformed using 
the Perdew-Wang LDA parametrization.t2l The TB cal- 
culations reproduce the LAPW results very well and are 
in good agreement with the experimental data. 

Phonons at T—0 K 

We determined the phonon dispersion curves and 
density-of-states (DOS) of fee Au by computing the dy- 
namical matrix. This was achieved by using a large su- 
percell, in our case containing 1331 atoms, and calculat- 
ing the forces on all atoms in response to the displace- 
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FIG. 1. Equation of State for selected crystal structures 
of gold using the tight-binding parameters discussed in the 
text. All coordinates are relaxed at each volume. The points 
are the first-principles LAPW energies used in the fit. From 
bottom to top, the ordering of the structures is fee (LAPW 
symbol 4-), hep, bee (LAPW symbol x), hexagonal co, A15, 
and simple cubic (sc, LAPW symbol ^(t). 

TABLE L Bulk modulus and elastic constants (in GPa) for 
Au computed using our TB model compared to the results of 
LAPW calculations and experimental data. All calculations 
are performed at the experimental room-temperature volume, 
the measured elastic copstants are taken from the compilation 
of Simmons and WangE3 
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ment of the atom at the origin. Provided this displace- 
ment is small enough it is possible to construct the real- 
space dynamical-matrix using finite differenceS|^nd com- 
pute the dynamical-matrix by a Fourier series.Q 

The high-symmetry direction phonon dispersion curves 
for Au at T=0 K-are shown in Fig. ^(a) together with ex- 
perimental data.E2l The overall structure of the dispersion 
curves is well reproduced. The low frequency transverse 
modes are in excellent agreement with experiment. The 
longitudinal higher-frequency modes, however, are sys- 
tematically too high close to the BZ edge. 

The phonon DOS is presented in Fig. ^(b) together 
with experimental results. The DOS has two main peaks. 
From to 3.5 THz the theoretical DOS reproduces the 
experimental data very well. In contrast in the region 
4-5 THz the position of the high frequency peak in the 
theoretical curve is overestimated by about 0.5 THz com- 
pared to experiment, consistent with the discrepancy in 
the frequency of the high-frequency longitudinal modes 
noted above. 

We checked that the observed discrepancies are not 
an artifact of the supercell method used to calculate 
the dynamical-matrix. In particular we made sure that 




FIG. 2. (a) Phonon-dispersion curves for Au at T=0 K, 
plotted along high-symmetry directions in the BZ. Lines are 
spline fits to the thaeretical TB data, open squares are experi- 
mental data points.cil (b) Phonon density-of-states (DOS) for 
Au at T=0 K. The full line is the theory result using our 
TB modeU-iwhile the dotted line is a fit to the experimen- 
tal values.B The dispersion curves and DOS were calculated 
from the dynamical-matrix computed using an fee supercell 
containing 1331 atoms. 

TABLE II. Selected phonon frequencies (in THz) at high 
symmetry points in the Brillouin-zone. The calculated TB 
and LAPW frequencies were obtained with the frozen phonon 
method in cells with 2 or 4 atoms. 
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an 11x11x11 supercell is large enough to converge the 
dynamical-matrix. A first check is given by the slopes of 
the dispersion curves as which are related to the 

elastic constants. We found that the slopes were consis- 
tent with our computed elastic constants (see Table |). 
Another check is to compute the frequency of BZ-edge 
phonons using the frozcn-phonon methocS for a 2 or 4 
atom unit cell, using the static code to calculate the 
total energies. The computed BZ-edge phonon frequen- 
cies are shown in Table ||. They are in perfect agree- 
ment with the frequencies derived from the dispersion 
curves obtained using the dynamical-matrix method (see 
Fig. ^) , hence confirming the accuracy of the latter ap- 
proach. Finally using the first-principles LAPW method 
we calculated the phonon frequencies at X and L and 
found very good agreement with experiment as shown in 
Table |l[ We therefore conclude that the overestimate of 
our calculated longitudinal phonon frequencies near the 
BZ-edge is a shortcoming of our TB parameters. An ob- 
vious approach to overcome this problem and improve 
the agreement with experiment would be to include the 
LAPW frequencies at X and L in our fitting database. 
Nevertheless we note that our TB results for the disper- 
sion curves present a substantial improvement over re- 
sults 0|btained using the second-moment approximation 
to TBE 

Electronic Density of States at finite temperature 
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FIG. 3. Temperature dependence of the electronic density 
of states of gold, calculated using the eigenvalues generated 
by the TBMD code, and averaged over ten time steps, as 
outlined in the text. 




FIG. 5. Finite-temperature phonon dispersion curves along 
the high-symmetry directions in the BZ, for Au at 300 K 
calculated using molecular-dynamics simulations based on 
our TB model. Filled circles are theoretical data, lines are 
polynomial fits to theoretical data, open squares are experi- 
mental data points.Ej The dispersion curves were computed 
by Fourier-transform of the time dependent wave-dependent 
VACF (see text). 
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FIG. 4. (a) Velocity-velocity auto-correlation functions 
(VACF) of Au calculated from molecular-dynamics simu- 
lations at 300 K and 1200 K using our TB model, (b) 
Finite-temperature phonon spectral-density (PSD) Z{f) for 
Au at 300 K and 1200 K. The PSD was computed by 
Fourier-transform of the VACF shown in (a) . 



The TBMD code computes all of the eigenvalues of the 
system at each time step, making it simple to determine 
electronic properties such as the density of states as a 
function of temperature. In Fig. ^ we show the electronic 
DOS at several temperatures. For each temperature, we 
saved the eigenvalues for ten different time steps. The 
DOS was then calculated assuming a Fermi distribution, 
and the resulting DOS were averaged. We see that the 
dominant effect of increasing temperature is to reduce 
the peaks in the electronic DOS spectrum. 

Phonons at finite temperature 

We determined the phonon dispersion curves and 
spectral-density of fee Au at finite temperature by per- 
forming MD simulations. In Fig. ^(a) we show the 
velocity- velocity auto-correlation functions (VACF) ob- 
tained from the MD simulation at 300 K and 1200 K. 
We see that increasing temperature damps out the oscil- 
lations in the VACF. 



The finite-temperature phonon spectral-density (PSD) 
can be— .obtained from the Fourier-transform of the 
VACF,E3| as shown in Fig. |(b). The PSD has the two 
well defined peaks, consistent with the DOS at T=0 K. 
The position of the peaks is also in agreement with the 
theoretical data at T=0 K. The limited resolution in the 
finite-temperature PSD, due to the short length of the 
MD simulation, makes a detailed comparison with ex- 
periment difficult. It should also be pointed out that 
the PSD is proportional to the phonon DOS only in an 
harmonic solid; one should therefore be cautious when 
comparing the PSD to the phonon DOS, in particular at 
high temperature. We believe this may explain the dif- 
ference in height between the two peaks in the PSD, in 
contrast to the T=0 K phonon DOS, where both peaks 
have about the same height. Comparison of the PSD at 
300 and 1200 K clearly reveals the effect of temperature: 
a clear shift of phonon frequencies to lower values and a 
broadening of the peaks in the PSD. 

To compute the phonon dispersion curves we calculate 
the PSD from the Fourier-transform of the velocity- and 
position-dependent auto-correlation function. Details p£ 
this computational procedure can be found elsewhere.cil 
The dispersion curves along high-symmetry directions in 
the BZ at 300 K are reproduced in Fig. |[ The com- 
parison with experimental data reveals the same discrep- 
ancies as in the T=0 K dispersion curves: the high fre- 
quency longitudinal modes are overestimated close to the 
BZ-edge. 

We determined the temperature dependence of phonon 
frequencies by performing MD simulations at 300, 600, 
900 and 1200 K where for each temperature we fix the 
volume at the experimental value. In Table III we show 



the frequency of selected BZ-edge phonons as a function 
of temperature. These frequencies were calculated as de- 
scribed above for the 300 K case. 

As expected the frequency of phonons decreases with 
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TABLE III. Selected phonon frequencies (in THz) at high 
symmetry points in the Brillouin-zone, calculated as a func- 
tion of temperature (T in Kelvin) from MD simulations using 
our TB model. 
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temperature. For the modes we computed the change 
was of the order of 0.2-0.3 THz, the transverse mode at 
L exhibiting the largest decrease. Part of this variation 
is probably due to increase in volume as temperature 
increases. It should be noted that using MD to com- 
pute the phonon dispersion curves at finite temperature 
can be particularly useful in systems where a particu- 
lar crystal structure is unstable at T=0 K [bcc Ti is one 
example) and where the dynamical-matrix method will 
predict unstable phonon modes. MD simulation may also 
be useful to determine the vibrational properties of sys- 
tems for which an experimental study may be difficult, 
e.g. clusters or nano-crystals. 

Thermal expansion 

To determine the theoretical thermal expansion coeffi- 
cient a we use the following definition for a: 
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This definition requires the calculation of the pressure as 
a function of temperature for a fixed volume. We perform 
MD simulations at 300, 600, 900 and 1200 K, keeping the 
volume fixed at the experimental value at room temper- 
ature (lattice constant a — 4.08 A). For each tempera- 
ture we selected 10 independent configurations from the 
trajectories generated by the MD and computed the in- 
stantaneous pressure. We found that 10 configurations 
per temperature were enough to get the average pressure 
with an error margin of ~5%. If, in Eq. |^, we assume 
the pressure varies linearly as a function of temperature 
and if for B we use the theoretical value of the bulk mod- 
ulus at T=0 K and at the experimental volume, we get 
a = 11 X 10~^ K~^. This underestirpates the experimen- 
tal value of 14x10-6 K"^ at 300 KEa 
An alternative definition of a is given by: 
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dV 
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To check the calculation of a based on Eq. |^ we computed 
a using this latter definition. This requires MD simu- 
lations for several volumes (typically 3 or 4) for a given 
temperature. For each volume V we compute the average 
pressure P. The equilibrium volume at each temperature 
is found by interpolating P{V) to find the volume that 
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FIG. 6. Lattice constant of An as a function of tem- 
perature. The black circles are the results of the molecu- 
lar-dynamics simulations lising our TB model, open squares 
are the experimental data.a 
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FIG. 7. Mean-square displacement of Au as a function of 
temperature. The filled circles are the results of the molecu- 
lar-dynamics simulations usmg our TB model, empty squares 
are the experimental pointsEZI 



gives zero pressure. We performed this procedure at 300, 
600, 900 and 1200 K to find V{T). 

In Fig. H we show the lattice constant as a function of 
temperature as deriveciiom the simulations, compared 
to experimental results.E3 The overall agreement with ex- 
periment is good given that the theoretical data is well 
within 1% of experiment in the temperature range we 
simulated. 

From Fig. ^ we can see that it is reasonable to assume 
that the volume varies linearly as a function of T. So by 
using Eq. ^, we get a = 11 x 10"^ K"^, in agreement 
with our previous estimate based on Eq. nl 

Mean-square displacement 

We used the atomic positions generated by the MD 
simulations performed for several temperatures at the 
corresponding experimental lattice constants to compute 
the atomic mean-square displacement (MSD). In Fig. 
we compare the temperature dependence of our com- 
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8. Pair correlation function g{r) of liquid gold at 



1773 K. The dotted line is the result of molecular-dynamics 
simulations usitig our TB model, the full line is the experi- 
mental result.cJ 

puted MSD with experimental data.0 The agreement 
with experiment is excellent up to 900K, at higher tem- 
peratures the theoretical MSD gets larger than the ex- 
perimental values. Again it should be noted that the 
results of our calculation are in much better agreement 
with experiment compared to previous work using the 
second-moment-approximation to TB.Id 

Liquid Au 

To further test the transferability of our model to dif- 
ferent atomic environments we studied the liquid phase 
of Au. The simulation was performed with 200 atoms in 
a periodic fee unit cell. The density of the sample was 
chosen to be eqrial to the experimental value of 16.746 g 
cm^-^ at 1773 K.Ea One thousand (1,000) MD steps were 
used to equilibrate the system. Statistical averages of 
structural properties were computed from data collected 
from the next two thousand (2,000) MD steps. The radial 
distribution function g(r) obtained from our simulation 
(shown in Fig. |^) is found to be in very good agreement 
with experimental dataHj 

We also calculated the electronic DOS of liquid gold at 
1773 K, using the same procedure as in the sohd. The 
result is shown in Fig. ^ While the overall shape of the 
DOS, including the width, is similar to Fig. ^, the high 
temperature and loss of symmetry has destroyed most of 
the peak structure. However, two new peaks appear at 
low energies, probably due to the lack of periodicity in 
the liquid. 



IV. CONCLUSIONS 

We presented results of simulations of bulk fee Au 
using our tight-binding model. Our TB Hamiltonian 
was used to compute the elastic constants of bulk Au, 




Energy (eV) 

FIG. 9. Electronic density of states of liciuid gold at 
T = 1773K, using the same method as in Fig. 



which were in very good agreement with the results of 
LDA calculations and experimental data. Using a super- 
cell method to compute the dynamical-matrix we deter- 
mined the phonon-dispersion curves and phonon density- 
of-states of Au at T=0 K. Our calculated dispersion 
curves are in good agreement with experimental data, 
except for a tendency to overestimate the frequency of 
longitudinal modes close to the Brillouin-Zone edge. We 
checked that this discrepancy is not a consequence of the 
method used to calculate the phonon frequencies. In ad- 
dition, we performed molecular-dynamics simulations at 
various temperatures to compute the phonon density-of- 
states and phonon-dispersion curves at finite tempera- 
ture. The molecular-dynamics simulation were also used 
to obtain the temperature dependence of the lattice con- 
stant and of the atomic mean-square-displacement. Both 
quantities were found to be in good agreement with ex- 
perimental data. Finally, we performed an MD simula- 
tion of the liquid phase of Au and obtained a radial distri- 
bution function in very good agreement with experiment. 
We believe these results demonstrate that our TB model, 
using parameters generated by the same procedurej^j can 
successfully be applied to the study of dynamical and fi- 
nite temperature properties of other metals. 
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